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INTRODUCTION 

The  LR  algorithm  is  a  method  for  determining  the  eigenvalues 
of  a  matrix.   The  basic  operation  utilized  by  the  algorithm  is 
the  triangularization  of  a  matrix  into  th   product  of  a  unit 
lower-triangular  matrix  and  an  upper-trie   ular  matrix.   A  unit 
lower-triangular  matrix  is  a  .atrix  whose  diagonal  elements  are 
all  1,  and  whose  elements  above  the  diagonal  are  all  zero.   An 
upper-triangular  matrix  is  one  whose  elements  below  the  diagonal 
are  all  zero. 

The  algorithm  will  first  be  introduced  in  its  simplest  form. 
Following  this  a  series  of  modifications  will  be  developed  to 
improve  the  accuracy  of  the  algorithm  and  to  accelerate  its 
convergence.   Gaussian  elimination  with  pivoting  is  an  elementary 
tool  in  developing  these  modifications.   This  process  is  a  basic 
tool  in  numerical  analysis  problems  dealing  with  systems  of 
equations  in  matrix  form.   A  development  of  Gaussian  elimination 
may  be  found  in  any  of  the  books  listed  in  the  bibliography. 


THE  LR  ALGORITHM 

The  algorithm  is  based  upon  the  triangular  decomposition  of 
a  matrix  A,  given  by 

[1]  A  s  LR 

where  L  is  a  unit  lower-triangular  matrix  and  R  is  upper- 
triangular.   If  we  now  form  the  similarity  transformation  L  AL 
on  the  matrix  A,  we  have 


[2]  L  1AL  =  L  X(LR)L  =  RL  . 


Hence  ,  if  we  decompose  A  and  then  multiply  the  factors  in 
reverse  order,  we  obtain  a  matrix  similar  to  A.   If  we  name  the 
original  matrix  A, ,  then  the  algorithm  is  defined  by  the 
equations 

™  As-1  =  Ls-lRs-l  '      Rs-lLs-l  =  As 

Thus  A   is  similar  to  A   -  and  by  induction,  to  A1 .   This 

process  is  repeated  until  we  obtain  a  matrix  Ag  such  that  Lg  =  I. 

which  means  the  diagonal  elements  of  Rg  are  the  eigenvalues  of 

A  .   Since  A  is  similar  to  A, ,  these  diagonal  elements  are  also 
s  s  1 

the  eigenvalues  of  A, .   This  then  is  the  LR  algorithm. 

Since  the  algorithm  is  based  upon  the  triangular  decomposi- 
tion of  a  matrix  A,  we  shall  introduce  a  method  for  the  triangu- 
larization  of  a  matrix.   For  the  original  matrix  A,  by  [1], 
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Theorem  1.   If  the  matrix  A  is  triangularized  such  that 
A  =  LR,  where  L  is  unit  lower-triangular  and  R  is  upper-triangular, 


:i 


det(A. . ) 
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ID 


det(A. . ) 
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,  i  =  1, 
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where  A. .  is  the  leading  principal  minor  of  A  of  dimension 


•  th 


(i)x(i),  and  A.,  denotes  this  minor  with  its  i    row  replaced  by 


•  th 
its  j    row 


Ji 


th 


If  we  partition  the  matrices  shown  above  along  the  i    row 

j  •  "th    , 
and  i    column,  we  have 
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It  follows 

that 

[4] 

A.  .  = 

L. -R. •  • 
11  11 

Let  A.,  denote  the  leading  (i)x(i)  principal 
3^ 

subma 

trix 

of 

A  with 

.  .   .th 
its  i   row 

,  ,   .,   .th 
'  replaced  by  its  3 

row .   Let  the 

same 

def 

ini 

tion 

hold  for  L. 

It  follows  that 
1 

A.  •  = 

31 

L. .R. •  • 
31  11 

Since 

L..  is  triangular,  so  is  L . . ,  but 

11                       31 

with 

k.  . 

31 

on 

the 

diagonal. 

Thus 

det(L. 

.)  =  k. - 
1     3i 

as  all  other  diagonal  elements 

are  1 .   Hence . 

1 

detCA^)  = 

k. -  detCR. • ) 
31      H 

. 

When  i  =  j 

detCA.^0 

=  det(R. • )  • 

Hence , 

det(A..) 
k..  =       31   , 

i  =  1 ,  2 , 

.  ,  n . 

31 

detCA. .) 

Similarly, 

we  can  find  an  expression  for  the 

r .  .  , 

13 

using 

Transposes 

t      1   t 

A. -  =  R. ;L. ■  , 

11    11  11 

t      t 
A .  .  =  R  •  •  L 

]i     31 

i 
ii 

i 

where  A   i 

s  the  transpose  of  A, 

Note  that 

Therefore 


detd!  .  )  =  1  . 


det(A..)  =  detCR.-)  ,     detCA..)  =  det(R!.) 

11         11   '         31         31 


T 

By  the  definition  of  R. . , 


and 


t 
,      r. .  detCR. . ) 
det(R..)  =  -2J ii- 

rii 


det ( R . . ) 
ii 


Therefore, 


so  that 


det(R.  .  .  ■_•) 
1-1,1-1 


r. .  detCR..)  det (R. . )det (R. . ) 

r..  =  -ii J1   =  ii 21— 

13      det(R.'.)  detCR.  .  .  ..)det(R!.) 

11  i-l,i-l      11 


det(A. . )det(A! . )  detCA^.) 

...  -         11 3i  -  31 

1D    det  (A.  ,  .  ,  )det(A.' .)  detCA.  ,  .  ,) 

i-l,i-l      11  i-l,i-l 
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Since  A. .  =  A. . , 

31   1: 


det (A. . ) 
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detCAi-l,i-l> 


where  det(AQ0)  =  1.   Thus 


det(A.  .).              .  det  (A.  .  ) 
r   =  3J ,    k..  =  ii_ 

"J    det  (A.    .   )       ^   det  (A..) 


For  a  simple  3x3  matrix  A, 


A   = 


Rll  a12  a13 


a21  a22  a23 


a31  a32  a33 


L   = 
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det(A      ) 


det(A31) 
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det(A22) 


R 
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det(A22) 


With  this  notation  established,  we  return  to  the  algorithm. 
The  following  assumptions  concerning  the  matrix  A  will  now  be 
made.   We  assume  the  eigenvalues  are  real  and  of  different  abso- 
lute value  and  that  the  leading  principal  minors  of  X  and  Y  are 
non-zero,  where  X  is  the  matrix  of  right-hand  eigenvectors  of  A 


and  Y  =  X 


,-1 


Theorem  2.   If  A  =  LR,  under  the  restrictions  just  stated, 


L  ■*■   I  anc 
s 


[5] 


n 


R  -*■  A  ■*> 
s     s 


X 


n 


as  S  -*■  «> 


where  X  denotes  the  possible  non-zero  elements  of  an  upper- 
triangular  matrix. 

To  prove  this  result,  we  establish  relations  between 
successive  iterations  which  will  be  used  extensively.   By  [3], 

Ac  =  f  ^   .L   -  ■. 

s     s-1  s-1  s-1 


Repeated  application  of  this  result  yields 
[6] 


-1   -1 

A   =  L   nL 
s 


-1  -1 
L0  Ln  A.,LnL, 


s-l-s-2  •'•  u2    ul   "1^2  •'•  Ls-2Ls-l 


or 


[7] 


L,L„  ...  L   _A„  =  A..L,  L„  ...  L 


Jl"2 


s-1  s     1  1^2 


Js-1 


Now  define  matrices  T   and  U   by 

s       s 


[8] 


"  L-L2 


L    and   U   =  R  R   , 
s        s    s  s-1 


J. 


These  matrices  are  unit  lower-triangular  and  upper-triangular 
respectively.   Consider  the  product  T  U  . 

T   U      =    L, L„    ...    L      , (L   R    )R      ,     ...    R0Rn 
ss  12  s-1      s    s      8—1  2    1 

=    L-i  L„    .  .  .    L      ,  A   R      -i     ...    R_R- 
1    2  s-1    s    s-1  2    1 

=   A-,  Ln  L0    ...     (L      -,  R      -,)     ...    R0Rt 
112  s-1    s-1  l    1 

=   A.,  L,  L„    ...    L      n  (A      .  )R      „    ...    Rr,Rn 
112  s-2      s-1      s-2  I    1 

2 

—      A-  Li-.  Ll0        ...        V  L  r.i\  rs  )        ...       K,-..!^-, 

112  S- I    S- Z  I    1 

2 

112  S-3        S-2        S-3  /     1 

3 

**     A -.  Li -,  L ,-.      ...      \  Lj        ",-K        /-.  y      ...     ix,.*  in.-, 

112  s-3    s-3  2    1 


=    AS 
Al 


Hence,  repeated  application  of  [7]  yields 


[9]  T  U   =  A^ 

s  s     1 


s 
so  that  T  U   gives  the  triangular  decomposition  of  A, 


PROOF  OF  THE  CONVERGENCE  OF  A 


Equation  [7]  is  a  fundamental  tool  in  the  analysis.   We 
shall  prove  that  if  the  eigenvalues  A.  of  A,  satisfy  the  relation 


\x\    >  ]A2|  >  |X3|  > 


n 


then  in  general  the  results  of  [5]  are  true. 

Before  giving  a  formal  proof,  we  consider  an  example  of  a 
matrix  of  order  three.   Note,  that  if  X~  AX  =  diagCX,.  ),  v/here  the 
X.    are  the  eigenvalues  of  A,  then  X  is  a  matrix  of  righthand 
eigenvectors  of  A,  written  as  column  vectors.   This  comes  from 
the  relation  AX  =  A.X,  where  X  denotes  the  eigenvector  corres- 
ponding to  X..   For  simplicity,  we  denote  the  matrix  X  of  right- 
hand  eigenvectors  in  the  form 


X  = 


x 


11 


:■: 


21 


>: 


12 


x 


22 


x31     X32 


"13 


x 


33 


and  its  inverse  Y  by 


X  -1  =  Y 


*U    Y12     y13 


y21   y22   y23 


*31     y32 


33 


If  we  denote  T  U   =  A,  =  B ,  then 
s  s     1 


10 


'11 


'12 


C13 


B  =  A*  =  X 


X 


'21 


'22 


'31 


'32      33 


where 


Cn  =  xixiiyn  +  x2xi2y2i +  A3xi3y3i 


c21  =  ^x2lYll  +  A^x22y21  ♦  X^x23y31 


s  s  s 

C31  =  AlX31yll  +  A2X32y21  +  X3X33y31 


c12  =  ^xllYl2  +  A^x12y22  +  A^x13y32 


C22  ~  AlX21y12  +  X2X22y22  +  X3X23y32 


'32 


XlX31y12  +  A2X32y22  +  X3X33y32 


'13 


XlXlly13  +  X2X12y23  +  X3X13y33 


'23 


XlX21y13  +  X2X22y23  +  X3X23y33 


C33  =  XlX31y13  +  X2X32y23  +  X3X33y33 


Now  T  U   is  the  triangular  decomposition  of  A-,  and  hence  the 

elements  of  the  first  column  of  T  will  correspond  to  the 

s  r 

elements  of  the  first  column  of  L  given  on  page  6.   They  are 


,    s       det(B   ) 

xll    x 


det  (b.,  ,  ) 


:  . 

,  (s) 
r21   " 

AlX21yll  +-  A2X22y21  +  X3X23^31  _  det^B2l) 
AlXllyll  +-  A2X12y21  +  A3X13y31    det^l^ 

,  (s) 

T31   " 

AlX31yll  +  A2X32y21  +  A3X33y31    ^^Sl5 

AlXllyll  +  X2X12y21  +  A3X13y31    det(Bn) 

Due  to  the  ordering  of  the  magnitude  of  the  eigenvalues,  \*   will 

dominate  the  denominator  as  s  ■*•  •».   Hence  if  x,  y    i    0, 

t21)  =  X21/Xll  +  °^2/Al)S 

t3?)  =  x31/xll  +  0.^2-/Xl')8 

where  0(A2/X  )s 

is  a  term(s)  of  order  (X  /X -}S,  which  approaches 

zero  since  | X,  j 

>  |x2|. 

Thus,  the  e 

lements  of  the  first  column  of  T   approach  the 

s 

corresponding  elements  in  the  triangular  decomposition  of  X. 

Similarly,  for 

the  elements  of  the  second  column  of  T  ,  we  have 

s  ' 

t(s)  =  1  __    det(B22) 

det(B22) 

[10]      t<3)  = 

det(B32)  _  (X1X2)Sa1  +  U^)^  +  ^2*3)% 

det(B22)    (X1X2)Sb1  +  (X1X3)Sb2  +  (X2X3)Sb3 

where 

al 

=  Cxnx32  -  5c31x12)(yily22  -  y21y12) 

a2 

=  (X11X33  "  X31X13)(ylly32  ~  y31y12) 

12 


a3    =    U12X33    "   X32X13)(y2ly32    "    y31y22} 
bl    =    (X11X22    "    X21X12)(ylly22    "    y21y12} 

b2  =  (x11x23  -  x21x13)(yily32  -  y31y12) 

b3  =  (x12x23  -  x13x22)(y21y32  -  y31Y22)   • 

If  we  now  divide  numerator  and  denominator  by  (X-^X^ 

s  s 

,  .    an  +  (X0/X0)  a„  +  (X„/X,)  a0 

(s  )    _1 3   2    2 3   1    3 

32     b1  +  (X3/A2)Sb2  +  (A3/X1)Sb3   ' 

Hence 

t(s)  5  *11*32  -  *31*12  t   o{ia/^,. 
XliX22  "  X21X12 

provided  (xi:lx22  "  X21X12  "*  ^ylly 22  ~  y21y12^  ^  ° '   We  Can  See 

(s  ) 

from  [10]  that  the  limiting  value  of  t„„   is  equal  to  the  corres- 
ponding element  obtained  in  the  triangular  decomposition  of  X. 
We  have  thus  established  that  if 

X  =  TU 
then,  provided  x11y1]_  ^  °  and  ^xllx22  "  X21X12  ^ylly  22  "  y2Iy12^ 

i   o, 

T   ->  T  . 
s 

Thus,  we  have  shown  in  this  simple  example  that  the  matrix 

T   tends  to  the  unit  lower-triangular  matrix  obtained  from  the 
s  to 

triangular  decomposition  of  the  matrix  X  of  eigenvectors , 


provided  the  leading  minors  of  X  and  Y  are  non-zero.   Wow  we 
shall  prove  Theorem  2  in  general  for  a  matrix  with  distinct 
eigenvalues .   If  we  write 


T  U   =  A,  s  B 
s  s     1 


then  we  have 


(s) 


t : 
31 


detCB.  .  )/det(B. . ) 
31       11 


as  we  have  proved  earlier  for  s  =  1 
From  the  relation 


we  have  that 


B  ■*■  Aj  ?  X  diagU?)  X  1 


B.  . 


[11] 


Xll      X12 


x21      x22 


xi-l,l   xi-l52 


Xjl      XJ2 


>; 


In 


'2n 


x.  . 

l-l  ,n 


x 


:n  J 


xiyn    xiyi2 


s        s 
A2y21    X2y22 


i  s 

a  y  1 


A  y  0 
n^n2 


A3>Ii 


^  s 

A2y2i 


1 s 

X   y  . 


1 


By  the  theorem  of  corresponding  matrices  (chapter  1,  section  15) 
det(B..)  is  equal  to  the  sum  of  the  products  of  the  corresponding 


1.   See  Bibliography 
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i-rowed  minors  of  the  two  matrices  on  the  right  in  [11].   Hence 
we  may  write 

Y    x(^}       y(i)       CA   A      A   )S 


tji   "  v-1  v(i)      ^i^ 


V  xU)      yU}       (A   A     A   ) 

Z_  PiP2.--pZPiP2-.-Pi  Pi  p2...  ?i 


where  x  J         is  the  i-rowed  minor  consisting  of  rows  1,  2, 
PlP2...Pi 

. . . ,  i-1  and  j  and  columns  p, ,  p0 ,  . . . ,  p .  of  X  and  y 

1    /         1  ^1-^2  *  '   i 

is  the  i-rowed  minor  consisting  of  rows  p, ,  p„,  ...,  p.  and 
columns  1,  2,     .  ..,  i  of  Y. 

The  dominant  terms  in  the  numerator  and  denominator  are 

s 
those  associated  with  (A_A„  ...A.)   providing  the  corresponding 

coefficients  are  non-zero.   Thus,  the  dominant  term  in  the 

denominator  is 

det(X..)  det(Y. . ) (X, A0  ...  A.)S 

11        11    1  2       1 

where  X..  and  Y . .  are  the  leading  principal  submatrices  of 

j.1       11  • 

order  i.   If  detCX..)  det(Y..)  is  non-zero,  we  have 

11        11 

,    v    det(X..)  detCY..)    det(X..) 
tCs)  ^      11 ii_  .      31 

3±  detCX..)  det(Y..)    detCX..) 

11       11         11 

showing  that  T   +  T,  where  T  =  XlT1.   We  have  from  [8],  [6],  and 
[7]  that   ' 

-I  -1        -1    -1 

A   =  T   .A,  T   ,  ■*  T   A,T  =  UX   AXU 

s     s-1  1  s-1        1 

,      v     -1 

=  U  diag(A, )  U 


showing  that  the  limiting  A   is  upper-triangular  with  diagonal 


elements  X . 


From  the  relation  L 


T   nT   and  equation  [12]  it  can  be 
s-1  s 


proved  in  a  similar  and  quite  tedious  argument  that 


kf?}  =  0(X./X.)s 

i]        1-  3 


as   s 


From  this  we  deduce,  using  the  relation  A   =  L  R  ,  that 

a.  .  =  OCX.  /A.  )S   as   6 '+  »  . 

13       i   3 

Hence,  if  the  separation  of  some  of  the  eigenvalues  is  poor,  the 
convergence  may  be  quite  slow. 

It  should  be  noted  that  in  establishing  these  results,  the 
following  assumptions  have  been  made,  either  explicitly  or 
implicitly . 

i) . ■  The  eigenvalues  are  real  and  of  different  magnitude. 

ii) .   The  triangular  decomposition  exists  at  every  stage. 
It  is  relatively  simple  to  construct  matrices,  not  otherwise 
exceptional,  for  which  this  is  not  true.   For  example, 


A 


0 
-3 


This  matrix  has  eigenvalues  1  and  3  but  has  no  triangular 
decomposition.   However,  this  case  can  be  handled  through  the 
use  of  interchanges j  which  we  will  introduce  later. 

iii)  .   The  leading  principal  minox^s  of  X  and  Y  are  all 
non-zero . 


16 


POSITIVE  DEFINITE  HERMIT IAN  MATRICES 

When  A.^  is  a  positive  definite  Hermitian  matrix  we  can 
remove  the  restrictions  of  the  last  section.   Thus 

A1  =  X  diag(Ai)  XH  ,  where  XH  =  X-1 

and  where  X  is  unitary  and  the  X.  are  real  and  positive.   Hence, 
equation  [10]  becomes 

y   x(j)      g(i)       (XX  X      )s 

(S)  _Z_,   p1p2...pi  Plp2.  .  .Pi  Plp2...Pi; 

Lji    T        K        I  I2  <*   *       X   )S 

z_,   P1P2- . - p±    Px  P2. . .  p£ 

Now  consider  the  class  of  aggregates  p1p2...p.  for  which 

Xp  p  ...p.  ^  °-   Let  ^1^2 ■ * ' ^i  be  a  member  of  this  class  such 

that  X   XQ     Xq   has  a  value  greater  than  that  of  any  other 
Hl  "2  * '  *i 

member.   Thus  the  denominator  is  clearly  dominated  by  the  term(s) 
in  (X   X      X   )s. 

As  for  the  numerator,  we  know  from  our  definition  of  a  a     a 

*1  ~2    -i 

that  no  term  exceeds  the  magnitude  of  (X   X     X   )s,  but  thi 


q±    q2. . .  q 


is 
l 


term  may  have  a  zero  coefficient.   Hence,  t)^)    tends  to  a  limit 
which  may  be  zero. 

No  assumptions  are  made  about  the  principal  minors  of  X  and 
hence  we  cannot  assert  that  the  limiting  T   is  the  matrix  obtained 
in  the  triangular  decomposition  of  X.   Accordingly,  we  proceed 
by  writing 

T   ->  T 
s     °° 


.  7 


so  that 


-1       -1 

S       S-l  S       oo    oo 


Further. 


[13]  A   =  T   nA,T   -  +  T   AnT 

S       S-l  1  S-l       oo    1  °° 


so  that  A   tends  to  a  limit,  say  A  .   Now 
s  °° 


-1     -1 

R   s  L   A  -*■    IT   AnT  , 

S       S    S         oo    1  oo  ' 


and  hence  R  tends  to  the  same  limit  as  A  .   Since  R   is  trian- 
s  s  s 

gular  for  all  s,  this  limit  must  be  triangular  also.   It  must 
have  the  eigenvalues  of  A,  in  some  order  on  its  diagonal,  since 
it  is  similar,  from  [13],  to  A,. 

Note  that  the  proof  is  unaffected  by  the  presence  of  multiple 
eigenvalues  or  by  the  vanishing  of  some  of  the  leading  principal 
minors  of  X,  though  the  A.  may  not  appear  in  decreasing  order  on 
the  diagonal  of  A  . 

It  is  now  advisable  to  assess  the  value  of  the  LR  algorithm 
as  a  practical  technique.   It  does  not  appear  to  be  very  adequate 
for  the  following  reasons : 

i).   Matrices  exist  which  have  no  triangular  decomposition, 
in  spire  of  the  fact  that  their  eigenproblem  is  well-conditioned. 
Without  some  sort  of  modification,  such  matrices  cannot  be 
handled  by  the  LR  algorithm.   Further,  there  is  a  much  larger 
class  of  matrices  whose  Triangular  decomposition  is  numerically 
unstable.   This  instability  can  arise  at  any  step  of  the  algorithm. 
which  may  lead  to  considerable  inaccuracy  in  the  computation  of 


lb 


the  eigenvalues.   A  method  to  eliminate  this  instability  will  be 
discussed  in  the  section  covering  inter-changes. 

ii).   The  volume  of  computation  is  very  high.   The  method 

2 
involves  (n-l)n  multiplications  at  each  step  of  the  process.   A 

method  of  reducing  the  amount  of  computation  will  be  discussed 

in  the  section  on  Hessenberg  matrices. 

iii).   The  convergence  of  the  subdiagonal  elements  of  L 

s 

depends  upon  the  ratio  C^r+1/Ap)  and  will  be  very  slow  if  separa- 
tion of  the  eigenvalues  is  poor.   This  problem  is  discussed  in 
the  section  covering  acceleration  of  convergence. 

Thus  if  the  LR  algorithm  is  to  be  useful,  it  must  be  modified 
to  meet  these  criticisms.   To  accomplish  this,  we  will  use 
elementary  matrices  which  interchange  or  combine  multiples  of 
the  rows  and  columns  of  A. 
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INTRODUCTION  OF  INTERCHANGES 

In  triangular  decomposition  numerical  stability  is 

maintained 

by  the  introduction  of  interchanges  if  necessary.   Consider  an 

analogous  modification  of  the  LR  algorithm. 

If  A  is  any  matrix,  by  Gaussian  elimination,  there 

exists  a 

product  of  elementary  matrices,  P,  such  that 

[14]                         PA  =  R 

where  R  is  upper-triangular.   We  can  therefore  complete 

a 

similarity  transformation  on  A  by  post-multiplying  R  by 

P"1  and 

have 

[15]                      PAP-1  =  RP"1  . 

In  particular,  when  there  are  no  interchanges  necessary 

the  matrix 

P  equals  the  matrix  L    and  [14-]  becomes 

L-1A  =  R  . 

In  this  case  the  matrix  on  the  right  of  [15]  is  RL. 

As  a  numerical  example  of  the  modified  process ,  we 

take  a 

3x3  matrix  to  illustrate  non-convergence  of  the  orthodox 

.  LR 

algorithm.   The  matrix  A,,  the  corresponding  A.,  and  the 

X  and  Y 

matrices  are: 

i      -l        i 

A,  =  5 

Al  = 

4     6-1 

A2  =  2 

4     4     1 

A3  =  1 

20 


X 


0    -1    -1 


1     0 


r 


l 


Y  = 


1     1 


-1 


z1      -1        h 


The  leading  principal  minor  of  X  is  zero.   Hence  we  cannot  be 
sure  that  the  orthodox  process  converges,  or,  if  it  does,  whether 
it  will  yield  the  eigenvalues  in  descending  order.   In  fact,  in 

this  case  the  elements  of  T   diverge.   In  Table  1,  we  have  shown 

s       °  ' 

the  results  of  the  first  three  steps  of  the  LR  process ,  and  from 

the  form  of  A  divergence  is  obvious . 
s 


1 

20 

4 


r'2 

•0.2 

6 

0.8 


1 

-5 

1 


100 


TABLE  1 
A. 


1   -0.04 


0.16 


•25 


1   -0.003 


500 


0.032 


-125 


Below:   LR  with  interchanges 


2 
3.2 


-i 


-1.25    1    1.25' 


1     0.8     1 


A3  A4 

5.157   3.040   -1.00   5.032   3.008   -1.00 


-0.174   1.833   1.042 


0.167       0.160      1.000 


-0.03      1.968      l.OOi 


0.032       0.032      1.000 


Ai 


A 


5 

3 

-1 

0 

2 

1 

0 

0 

1 

Tabie  1  also  gives  the  results  obtained  using  the  modified 
process.   The  first  step  in  detail  yields: 


1-11 
4     6-1 


R(l,2) 


1    -1 


-1 


R2-(1/4)R1 


4      6-1 


0    -2.5    1.2! 


where  R(l,2)  indicates  an  interchange  of  rows  1  and  2,  and 
R2-(l/4)Ri  indicates  the  subtraction  of  1/4  the  elements  of  row  1 
from  the  corresponding  elements  in  row  2 . 


-1 


0    -2.5    1.25 


R3-R1 


-2 


-1 


2.5    1.25 


R3-.8R2 


-1 


0    -2.5    1.25 


=  R 


-1 


"1 


0    -2.5    1.25 


CC1.2) 


-2.5 


1.25 


r 


C1+(1/4)C2 


-2.5     0     1.25 


< v 


C1  +  C3 


r 


-i 


-1.25    0     1.25 


C2+.8C3 


-1.25 


3.2 


-1 


1.25 


=  A, 


The  choice  of  which  elements  to  interchange  in  the  first 
step  is  an  important  one.   Since  we  have  two  equal  elements  in 
the  first  column,  the  second  row  is  interchanged  with  the  first 
This  choice  is  made  so  that  the  largest  diagonal  element,  6, 
will  become  the  leading  element  of  the  matrix  after  the  corres- 
ponding column  interchange.   This  is  done  because  we  want  the 
eigenvalues  in  descending  order  on  the  diagonal. 

Only  one  interchange  is  necessary  in  the  first  step.   In 
the  subsequent  steps  we  are  using  the  orthodox  LR  technique,  as 
no  more  interchanges  are  necessary.   The  matrix  A   converges 
very  rapidly  to  an  upper-triangular  matrix.   The  use  of  inter- 
changes has  not  only  yielded  convergence,  but  also  gave  the 
eigenvalues  in  descending  order. 


2  3 


The  use  of  interchanges  has  given  us  numerical  stability 
in  triangular  decomposition  for  a  simple  reason.   As  may  be 
noted  in  the  simple  example  on  page  6,  triangularization  utilizes 
a  division  by  det(A  .).   The  introduction  of  interchanges  has 
made  detCA,.,)  greater  than  or  equal  to  all  other  elements  in  the 
first  column,  hence,  the  first  and  second  columns  of  L  will 
become  increasingly  smaller.   Without  this  modification,  the 
orthodox  LR  technique  performs  this  division  without  regard  to 
the  relative  size  of  a,,,  and  hence,  as  is  seen  in  Table  1, 
elements  may  increase  in  size. 


THE  UPPER  HESSENBERG  FORiM 

It  would  seem  that  the  volume  of  work  involved  is  still 
prohibitive  when  A  is  a  matrix  with  few  zero  elements.   If, 
however,  A  is  of  a  condensed  form  which  is  invariant  with  respect 
to  the  LR  algorithm,  then  the  volume  of  work  might  well  be 
reduced. 

The  major  form  which  meets  this  requirement  is  the  upper 
Hessenberg  form,  which  is  invariant  with  respect  to  thex modified 
LR  algorithm,  and  therefore  a  fortiori  with  respect  to  the 
orthodox  process.   We  first  reduce  a  matrix  to  upper  Hessenberg 
form.   We  define  a  matrix  to  be  upper  Kessenberg  form  if  a..  =  0 
for  i  >_  j+2.   This  reduction  may  be  accomplished  in  the  following 
manner.   We  will  utilize  a  matrix,  call  it  N,  with  the  following 
characteristics:   N  is  unit  lower-triangular,  and  the  elements 
of  the  first  column,  except  for  the  diagonal  element,  are  zero. 
Hence,  N  has  the  form,  for  n  =  5, 


N  = 


0 


0     0 


0 

n32 

1 
J. 

0 

0 

0 

n42 

n43 

1 

0 

0 

n52 

n53 

n54 

1 

For  our  initial  matrix  A,  we  now  form  the  equation 


[16] 


AN  =  NH 


where  H  is  upper  Hessenberg 
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LAI 

X     X     X  X  X 

X     X     X  X  X 

X     X     X  X  X 

X     X     X  X  X 

X     X     X  X  X 


1     0 


0    n 


32 


EN] 
0 

0 


0    n42   n43 


0    n52   n53   n54 


0 

0 

0 

0 

0 

0 

1 

0 

54 

1 

[N] 

0     0     0 

10     0 


0    n 


32 


1     0 
1 


0    n42   n43 


°    n52   n53   n54 


[H] 

0   h   h   h  h  h 

0   h   h   h  h  h 

0   0   h   h  h  h 

0   0    0   h  h  h 

10    0    0  h  h 


where  x  and  h  denote  possible  non-zero  elements.   The  elements 
of  N  and  H  can  be  computed  column  by  column.   Multiplication  of 
[16]  by  N   yields 


N  1AN  =  H 


Hence,  A  is  similar  to  H,  an  upper  Hessenberg  matrix.   For  a 

further  discussion  of  this  method,  see  Wilkinson  (chapter  6, 

2 

section  11)  . 


2.   See  Bibliography 
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We  now  return  to  the  invariance  of  the  upper  Hessenberg 
form  with  respect  to  the  LR  algorithm.   We  want  to  show  that 
after  one  complete  step  of  the  algorithm  A?  is  of  upper  Hessen- 
berg form.   This  is  done  by  induction.   This  can  be  done  in  this 
manner  since  the  elementary  operations  used  to  triangularize  A, 
will  involve  only  rows  1  and  2  at  step  1,  only  rows  2  and  3  at 
step  2,  and  so  forth.   Hence,  after  (r-1)  steps  of  the  post- 
multiplication  by  these  factors,  only  the  first  (r-1)  columns 
of  the  matrix  will  have  been  affected.   Thus,  assume  that  after 
(r-1)  steps  the  matrix  is  of  upper  Hessenberg  form  in  its  first 
(r-1)  columns  and  triangular  in  the  remaining  columns.   For  the 
case  n  =  6 ,  r  =  4 ,  it  is  of  the  form 


X     X     X     X     X     X 


X     X     X     X     X     X 


X     X     X     X     X 


X     X     X     X 


X     X 


X 


where  x  denotes  possible  non-zero  elements,  and  the  elements 
not  shown  are  ail  zero.   The  next  step  will  be  the  interchange 
of  columns  r  and  (r+1)  if  rows  r  and  (r  +  1)  were  interchanged;  if 
not,  there  is  no  effect.   The  resulting  matrix  is  therefore  of 
the  form  (a)  or  (b)  shown  below. 


(a) 


X     X     X     X 


X     X 


X 


X     X     X     X 


X 


X     X     X     X 

x    0    x 

X 


(b) 


X     X 


X 


X    X 


X    X 
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X     X     X     X     X     X 


X     X     X     X 


X     X 


X 


where  (a)  shows  the  effect  of  an  interchange. 

Next  a  multiple  of  column  (r+1)  will  be  added  to  column  r 
and  the  resulting  matrix  will  be  of  form  (c) ,  from  (a),  or  (d) , 
from  (b) . 


(c) 


X     X     X     X     X 


X 


X     X     X     X  X  X 

X     X     X  X  X 

X     X  X  X 

x  0  x 


(d) 


X     X 


X     X     X     X 


X     X 


X     X     X     X     X 


X     X     X     X 


X     X 


X 

In  either  case  the  matrix  is  of  upper  Hessenberg  form  in  its 
first  r  columns  and  Triangular  otherwise;  the  extra  zero  element 
in  (c)  being  of  no  significance.   Hence,  the  upper-Hessenberg 

form  is  invariant  with  respect  to  the  modified  algorithm. 

2 
There  are  n  /2  multiplications  in  the  reduction  to  triangular 

form  and  n  /2  in  the  post-multiplication,  yielding  n2  steps  in 

one  complete  cycle  of  the  LR  algorithm,  compared  with  (n-l)n2 

for  a  full  matrix. 
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ACCELERATION  OF  CONVERGENCE 

Although  a  preliminary  : 

reduction  to  upper  Hessenberg  form 

reduces  the  volume  of  comput 

ation  considerably,  the  modified  LR 

method  will  still  be  somewha 

t  uneconomical  without  improving 

the  rate  of  convergence. 

As  we  have  seen  for  a  g> 

sneral  matrix,  the  elements  in 

positions  (i,j),  i  >  j  ,  tend 

s 
to  zero  as  (A. /A.)   does.   For 

i  : 

Hessenberg  matrices  the  only 

non-zero  subdiagonal  elements  are 

those  in  positions  (i+l,i). 

Now  consider  the  matrix  (A  -  pi), 

where  p  is  some  real  number. 

We  will  discuss  a  method  of 

selecting  p  later. 

The  matrix  (A  -  pi)  has 

eigenvalues  (A.  -  p)  and,  for 

(s) 

example,  the  element  a     ,  ' 

n,n-l 

tends  to  zero  as  [(A   -  d)/(A   -  -  p] 

n        n-1 

does .   If  p  is  a  good  approx 

imation  to  A  ,  the  element 

n ' 

a     .,  will  cecrease  rapidly 
n,n-i                 *    J 

Thus  it  would  be  to  our  advantage 

to  use  (A   -  pi)  rather  than 
s 

A  . 
s 

Note,  in  particular,  if 

p  is  exactly  equal  to  A  ,  then  the 

(s) 

element  a     .  would  be  zero 
n  ,n-l 

after  one  iteration.   This  may  be 

seen  by  considering  the  tria: 

ngularization  process .   None  of  the 

pivots  (the  diagonal  element; 

3  of  the  matrix  R)  can  be  zero 

except  the  last,  because  at  i 

aach  stage  in  the  reduction  the  pivot 

is  either  a.  +,  -or  some  oth< 

5r  non-zero  number  and  we  are  assum- 

ing  that  no  a.+.  .  is  zero  originally.   Hence,  since  the  deter- 

minant of  (A  -  X  I)  is  zero, 

n           ' 

and  R  is  the  form 

X 


X 


X 


X 

X 

X 

X 

X 

X 

X 

X 

0_ 

the  whole  of  the  last  row  being  null.   The  subsequent  post- 
multiplication  merely  combines  columns  and  after  the  iteration 
the  matrix  is  of  the  form 


X     X     X     X 


X     X 


X 


X     X 


X 


XXX 

0    0 


The  previous  discussion  suggests  the  following  modification 

th 

of  the  LR  method.   At  the  s  J  stage  we  perform  a  triangulr 

decomposition  of  (A^  -  k  I) ,  where  k   is  some  suitable  value, 

s  s  s  ' 

rather   than   of   A    .      We   therefore   produce   the    sequence   of  matrices 

s  F  ■  ^ 

defined   by 


A      -   k    I    =    L   R 
s  s  s    s 


RL      +    k    I    =   A   ^, 
s    s  s  s+1 


Thus 


A    , 
s+ 


R   L      +   k   I    =    L   X(A      -   k    I)L      +   k   I    =    L   TA   L 
ss  s  ss  ss  s  sss 


and  hence  the  matrices  A  are  again  similar  to  A,.   In  fact 

\  =  L_1A  L   =  if1!/"1,  A   nL   , Lo  =  L~   . . .  L~  L~  A-.L-.Lj  ...  L 

s+1     s   3  s     s   s-1  s-1  s-1  S     S        S   1   X  J.  I 


ov 


L1L2  •*•  LsAs+l  =  A1L1L2  '••  Ls  ' 

This  formulation  of  the  modification  is  usually  described  as  LR 
with  shifts  of  origin  and  restoring,  because  the  shift  is  added 
back  at  each  stage.   We  have 

L1L2  •••  Ls-l(LsVRs-l  '••  R2R1 

-    L1L2  •••  Ls-l(As  "  ksI)Rs-l  •'•  R2R1 
=  (A1  -  ksI)LlL2  ...  (L^R^)  ...  Vi 

=  (A1  -  ksI)(A1  -  ks_1I)L1L2  ...  Ls_2Rs_2  ...  R^ 
=  (A,  -  k  I)(An  -  k  .1)  ...  (A,  -  k  I)  . 

1      S      1      S-1  1      -L 

Hence ,  writing 


Ts  =  LlL2  ...  Ls,      Us  =  Rs  ...  R^ 

we  see  that  T  U   gives  the  triangular  decomposition  of 
s  s 


n  (ai  -  kin 


i=l 

the  order  of  the  factors  being  immaterial. 

In  a  practical  sense  we  are  now  faced  with  selecting  a 
suitable  sequence  of  k   so  as  to  give  rapid  convergence.   We 


'.. 


(s)         (.<=■ ) 
expect ,  in  the  matrix  A  ,  that  a    .  and  a    will  approach 

s         n,n-l      nn        ri 

zero  and  A   respectively.   Hence  it  is  reasonable  to  take  k   = 

n    c  J  s 

a„„   as  soon  as  a     .  becomes  <  1  or  a     indicates  it  is 
nn  n,n-l  nn 

C  s  ) 

converging.   In  fact,  it  is  simple  to  show  that  when  a    ,  is 

^  n-n-1 

of  order  £  then,  if  we  choose  k   =  a    ,  we  have 

s    nn  ' 


(s+1)    n,    2, 
a    .  =  0(e  )  . 
n,n-l 


(s). 


For  (A  -  a^'i)  with  n  =  6,  is  of  the  form  shown  in  matrix  (a) 

O  ill  1 

below . 


(a) 


X  X  X  X 


X  X  X  X  X 


X  X 


X 


XXX 
XXX 
XXX 

£  0 


(b) 

X            X            X            X            X  X 

X            X            X            X  X 

X            X            X  X 

XXX 

a  b 

e  0 


(c) 

X  X  X  X  X 

X  X  X  X  X 

XXX  X 

XX  X 

a  b 
-be/c 


Consider  now  the  reduction  of  (A   -  a    I)  to  triangular  form 

s     nn 

by  the  use  of  Gaussian  elimination  with  interchanges.   Matrix  (b) 
above  indicates  the  form  of  the  matrix  when  only  the  last  row 
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remains  to  be  reduced.   The  element  (a)  in  row  (n-1)  will  not 

(s ) 

be  small  unless  a    happened  to  bear  some  special  relationship 

to  the  leading  principal  minor  of  order  (n-1)  of  A  ,  i.e.,  if 

s       ' 

the  shift  was  an  eigenvalue  of  this  minor.   Then  no  interchange 
would  be  necessary  in  the  last  step  and  we  would  take  row  (n) 
-  e/a  row  (n-1),  i.e.,  R(n)  -  e/a  R(n-l).   The  triangular  matrix 
is  thus  of  form  (c)  above. 

When  we  have  post  multiplied  by  all  factors  except  those 
involving  the  last  two  columns ,  the  current  matrix  will  be  of 
form  (a)  or  (b)  given  below,  depending  on  whether  an  interchange 
did  or  did  not  occur,  i.e.,  C(n-l,n). 


x 


(a) 


x        x 


XXX 


X 


X  X  X  X  X  X 


X  X  X  X 


X 


0  b 


■be/a. 


x        x 


(b) 


x 


x 


X  X  X  X  X 


X  X  X  X 


XXX 


x        a 


x 


X 


X 


X 


-be/a 


To  complete  the  post-multiplication,  we  add  e/a  times  column  n 
to  column  (n-1),  no  interchange  being  necessary  in  general. 
The  final  matrix  is  of  form  (a)  or  (b)  below. 


(a) 


(b) 


XXX 


X     X     X     X 


XXX 


X     X 


X 


X 


X 


X 

be/a 

-b£2 


X 


X 


X 


-bs 


X     X     X 


X 


X 


X     X 


X 


X 

a 
-be2 


X 


Hence,  axter  restoring  the  shift,  we  have 


(s+1)     (s)    ,   , 
a       =  a      be/a, 

nn       nn         ' 


(s  +  1) 

i    , 

n  ,n-l 


-  2,2 

-be  /a 


,  ^   (s+1)     .  ,  .      ,   (s+1)  .    -    .    2 

so  tndi.  a      is  lnaeea  converging  and  a     ,  is  or  order  e  , 
nn  °  &      n ,n-l  ' 

which  we  wished  to  show.   Note  that  any  interchange  which  may 
take  place  in  the  other  steps  of  the  reduction  are  of  little 

significance . 

(s) 

in  general,  once  a     1    has  become  small,  it  will  diminish 
°  n,n-l  ' 

rapidly  in  value.   When  it  is  negligible  to  working  accuracy  we 

(s) 

can  treat  it  as  zero  and  The  current  value  of  a     is  then  an 

nn 

eigenvalue.   The  remaining  eigenvalues  are  those  of  the  leading 

principal  submatrix  of  order  (n-1).   This  matrix  is  itself  of 

Hessenberg  form  so  that  we  can  continue  with  the  same  method, 

working  with  a  matrix  of  order  one  less  than  the  original. 

Since  we  are  expecting  all  sub-diagonal  elements  to  tend  to  zero, 
(s) 


n 


,    2  ^^y  already  be  fairly  small,  and  in  this  case  we  can 


immediately  use  a  ,    ,  as  the  next  value  of  k  . 
J  n-l,n-l  => 

Continuing  in  this  way  we  may  find  the  eigenvalues  one  by 
one,  working  with  matrices  of  progressively  decreasing  order. 
The  later  stages  of  the  convergence  to  each  eigenvalue  will  be 
quadratic  generally,  and  moreover,  we  can  expect  that  when 
finding  the  later  eigenvalues  in  the  sequence,  we  shall  have  a 
good  start. 


CONCLUSION 

The  basic  LR  algorithm,  although  it  introduces  the  basic 
method  of  operation,  has  been  seen  to  be  quite  tedious  and  quite 
possibly  inaccurate.  However,  the  introduction  of  interchanges, 
reduction  to  upper  Hessenberg  form,  and  shifts  of  origin  have 
given  numerical  stability,  a  reduction  of  the  volume  of  computa- 
tion, and  accelerated  convergence  respectively.  These  modifica- 
tions have  made  the  LR  algorithm  a  somewhat  practical  method  for 
the  computation  of  the  eigenvalues  of  a  matrix  with  the  restric- 
tions that  were  placed  upon  it. 

These  restrictions,  mainly  that  the  matrix  have  real 
distinct  eigenvalues,  are  still  quite  prohibitive.   It  is 
generally  quite  difficult,  if  not  impossible,  to  determine 
whether  the  matrix  A  has  distinct  real  eigenvalues.   There  exist, 
however,  more  advanced  methods  which  are  able  to  handle  matrices 
with  these  restrictions.   The  QR  algorithm,  which  is  analogous 
to  the  LR  technique,  is  more  adaptive  to  these  limitations. 

For  a  discussion  of  the  QR  algorithm,  see  Wilkinson  (chapter  8, 

3 
section  28 ) . 


3 .   See  Bibliography 
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The  LR  algorithm  is  an  iterative  method  for  determining  the 
eigenvalues  of  a  matrix.  The  basis  of  the  algorithm  is  the 
triangularization  of  a  matrix  A,  into  the  product  of  a  unit  lower- 
triangular  matrix  L  and  an  upper-triangular  matrix  R.  This  yields 
the  equation  A,  =  LR.  Multiplication  of  these  matrices  in  reverse 
order  yields  a  matrix  similar  to  A,,  i.e.,  RL  =  A2 .  The  algorithm 
is  thus  defined  by  the  equations 

A   .  =  L   .R   . 
s-1     s-1  s-1 

R   nL   .  =  A   . 
s-1  s-1    s 

As  s  approaches  infinity,  L  ■*■  I  and  R  ■»■  A  ...   Hence,  the  desired 
eigenvalues  of  A-  are  the  diagonal  elements  of  R  .   This  is  true 
in  general,  however,  only  if  the  eigenvalues  are  distinct  and 
real.   This  orthodox  procedure  can  be  numerically  unstable  and 
quite  slow  to  converge,  unless  some  modifications  are  made. 

The  introduction  of  Gaussian  elimination  with  interchanges 
eliminates  the  problem  of  numerical  instability  by  assuring  That 
leading  principal  submatrices  of  A  are  non-zero. 

The  process  still  requires  a  large  amount  of  computation. 
To  relieve  this  problem,  the  matrix  A  is  reduced  to  an  upper 
Hessenberg  matrix.   The  upper  Hessenberg  matrix  is  invariant  with 
respect  to  the  LR  algorithm.   This  reduces  the  computation  from 
(n-l)n   multiplications  to  nl   multiplications  per  step. 

The  problem  of  the  possible  slow  rate  of  convergence  of  the 
algorithm  is  remedied  by  the  use  of  shifts  of  origin.   By  working 
with  the  matrix  (A  -  pi),  where  p  is  an  approximation  to  A . ,  the 


rate  of  convergence  can  be  considerably  improved.   When  the  sub- 
diagonal  element  (a.,..)  is  less  than  1,  choose  p  =  a... 
These  modifications  definitely  improve  the  practical 
application  of  the  basic  algorithm,  however,  two  major  restric- 
tions were  assumed  in  the  development.   It  was  assumed  that  the 
matrix  A  had  real  and  distinct  eigenvalues.   These  two  limitation: 
are  difficult  to  recognize  in  most  problems ,  and  they  require 
additional  methods  to  handle  them.   However,  if  these  limitations 
are  not  present,. the  modified  LR  algorithm  will  yield  the  eigen- 
values of  the  matrix  A,  and  the  method  will  be  practical  to  use. 


